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SUMMARY 

New elements are derived to validate and assess the assumed natural deviatoric strain (ANDES) formulation. 
This is a brand new variant of the assumed natural strain (ANS) formulation of finite elements, which has 
recently attracted attention as an effective method for constructing high-performance elements for linear 
and nonlinear analysis. The ANDES formulation is based on an extended parametrized variational principle 
developed in recent publications. The key concept is that only the deviatoric part of the strains is assumed 
over the element whereas the mean strain part is discarded in favor of a constant stress assumption. Unlike 
conventional ANS elements, ANDES elements satisfy the individual element test (a stringent form of the 
patch test) a priori while retaining the favorable distortion-insensitivity properties of ANS elements. The 
first application of this new formulation is the development of several Kirchhoff plate bending triangular 
elements with the standard nine degrees of freedom. Linear curvature variations are sampled along the three 
sides with the corners as *gage reading* points. These sample values are interpolated over the triangle using 
three schemes. Two schemes merge back to conventional ANS elements, one being identical to the Discrete 
Kirchhoff Triangle (DKT), whereas the third one produces two new ANDES elements. Numerical experiments 
indicate that one of the ANDES element is relatively insensitive to distortion compared to previously derived 
high-performance plate-bending elements, while retaining accuracy for nondistorted elements. 
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1. INTRODUCTION 

Despite almost three decades of work, plates and shells remain a important area of research in finite element 
methods. Challenging topics include: 

1. The construction of high performance elements. 

2. The modeling of composite and stiffened wall constructions. 

3. The treatment of prestress, imperfections, nonlinear, dissipative and dynamic effects. 

4. The development of practical error estimators and adaptive discretization methods. 

5. The interaction with nonstructural components, for example external and internal fluids. 

This paper addresses primarily the first challenge, although it must be recognized that progress in this 
direction is shaped to some extent by thinking of the others. The main motivation here is the construction 
of simple and efficient finite elements for plates and shells that are lock-free, rank sufficient and distortion 
insensitive, yield accurate answers for coarse meshes, fit naturally into displacement-based programs, and 
can be easily extended to nonlinear and dynamic problems. Elements that possess these attributes to some 
noticeable degree are collectively known as high performance or HP elements. 

Over the past three decades investigators have resorted to many ingenious devices to construct HP elements. 
The most important ones are listed in Table 1. The underlying theme is that although the final product 

may look like a standard displacement model so as to fit naturally into existing finite element programs, 

the conventional displacement formulation is abandoned . (By “conventionar we mean the use of conforming 
displacement assumptions into the total potential energy principle.) 

1.1 A Unified Variational Framework 

Table 1 conveys the feeling of a bewildering array of tools. The question arises as to whether some of them 
are just facets of the same thing. Limited progress has been made in this regard. One notable advance in 
the 1970s has been the equivalence of reduced/selective integration and mixed methods achieved by Malkus 
and Hughes [25], 

The present work has benefited from the unplanned confluence of two unification efforts. An initial attempt 
to place the free formulation [6-9] within the framework of parametrized hybrid variational principles was 
successful [14-16]. The free formulation in turn “dragged” incompatible shape functions, the patch test, 
and energy balancing into the scene. Concurrently a separate effort was carried out to set out the assumed 
natural strain (ANS) and projection methods in a mixed/hybrid variational framework [26,27]. Comparison 
of the results led to the rather unexpected conclusion that a parametrized variational framework was able 
to encompass ANS and the free formulation as well as some hitherto untried methods [17,18]. 

The common theme emerging from this unification is that a wide class of HP elements can be constructed 
using two ingredients: 

(1) A parametrized functional that contains all variational principles of elasticity as special cases. 

(2) Additional assumptions (sometimes called ‘Variational crimes” or “tricks”) that can be placed on a 
variational setting through Lagrange multipliers. 

As of this writing it is not known whether the “wide class” referred to above encompasses all IIP elements or 
at least the most interesting ones. Some surprising coalescences, such as DKT and ANS bending elements, 
however, have emerged from this study. 
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Table 1* Tools for Constructing HP Elements 



Technique 

Year introduced 

1 . 

Incompatible shape functions 

early 1960s 

2. 

Patch test 

1965 

3. 

Mixed and hybrid variational principles 

1965 

4. 

Projectors 

1967 

5. 

Selective reduced integration 

1969 

6. 

Uniform reduced integration 

1970 

7. 

Partial strain assumptions 

1970 

8. 

Energy balancing 

1974 

' 9. 

Directional integration 

1978 

10. 

Limit differential equations- 

1982 

11. 

Free formulation 

1984 

12. 

Assumed natural strains 

1984 


1.2 The Assumed Natural Strain Formulation 

The assumed natural strain (ANS) formulation of finite elements is a relatively new development. A restricted 
form of the assumed strain method (not involving natural strains) was introduced in 1969 by Wiliam [33], who 
constructed a 4-node plane-stress element by assuming a constant shear strain independently of the direct 
strains and using a strain-displacement mixed variational principle. (The resulting element is identical to that 
derivable by selective one-point integration.) A different approach advocated by Ashwell [1] and coworkers 
viewed “strain elements” as a way to obtain appropriate displacement fields by integration of appropriately 
assumed compatible strain fields.* 

These and other forms of assumed-strain techniques were overshadowed in the 1970s by developments in 
reduced and selective integration methods. The assumed strain approach in natural coordinates, however, 
has recently attracted substantial attention [2,11,19,22,24,29,30,31], particularly in view of its effectiveness 
in geometrically nonlinear analysis. 

As noted above, the unification presented in [17,18] merges two HP element construction schemes: the 
free formulation (FF) of Bergan and Nygard [8] and a variant of ANS called ANDES (acronym for Assumed 
Natural Deviatoric Strains) described in further detail below. The stiffness equations produced by the unified 
formulation enjoy the fundamental decomposition property summarized in Box 1. 

In the ANDES variant of ANS, assumptions are made only on the deviatoric portion of the element strains, 
namely that portion that integrates to zero over each element. This assumption produces the higher order 
stiffness labeled K *22 in Box 1. The mean portion of the strains is left to be determined variationally from 
assumptions on the limit stress field, and has no effect on the stiffness equations. 


* In fact, this was the technique originally used by Turner et al . [32] for deriving the constant-strain membrane 
triangle in their celebrated 1956 paper. 
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This paper describes the construction of the first ANDES elements. These are Kirchhoff plate-bending 
triangular elements with the standard 9 degrees of freedom (one displacement and two rotations at each 
corner). This choice is made because of the following reasons: 

1. High-performance three-node triangular plate bending elements, whether based on Kirchhoff or Fteissner- 
Mindlin mathematical models, have not been previously obtained through the ANS formulation. (Al- 
though the DKT element [3] qualifies as high-performance and^ is in fact an ANS element as shown 
later, it has not been derived as such.) The situation is in sharp contrast to four-node quadrilateral 
bending elements, for which HP elements have been constructed through a greater variety of tools; see 
e.g. [11,20,23,24,29]. 

2. High performance elements of this type have been obtained through the FF and ancestors of the FF 
[6,7,8,12], and they are considered among the best performers available. It is therefore intriguing whether 
elements based on the ANDES variant can match or exceed this performance. 

The basic steps in the construction of K* and K* for a general three-dimensional element are summarized 
in Boxes 2 and 3, respectively. For justification of these “recipees” the reader is referred to [17,18]. The 
derivation of conventional ANS elements is summarized in Box 4. 
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Box 1 Decomposition of the Element Stiffness Equations 

Let K be the element stiffness matrix, v the visible element degrees of freedom (those 
degrees of freedom in common with other elements, also called the connectors) and f the 
corresponding element node forces. Then the element stiffness equations decompose as 

Kv = (K* + K A )v = f. (1) 

K& and K* are called the basic and higher order stiffness matrices, respectively. The basic 
stiffness matrix, which is usually rank deficient, is constructed for convergence. The higher 
order stiffness matrix is constructed for stability and (in more recent work) accuracy. A 
decomposition of this nature, which also holds at the assembly level, was first obtained by 
Bergan and Nygard in the derivation of the free formulation [8]. 

In the unified formulation presented in [17,18] the following key properties of the decompo- 
sition (1) are derived. 

1. Kj is formulation independent and is defined entirely by an assumed constant stress state 
working on element boundary displacements. No knowledge of the interior displacements 
is necessary (Box 2). The extension of this statement to C° plate and shell elements is 
not straightforward, however, and special considerations are necessary in order to obtain 
K& for those elements. 

2. K* has the general form 

K* = jttKhtt + J22Ka 22 + j23Kft23- (2) 

The three parameters j 2 2 , ji 3 and j 33 characterize the source variational principle in 
the following sense: 

(a) The FF is recovered if j 22 = j23 = 0 and J 33 = 1 — 7, where 7 is a K* scaling 
coefficient studied in [9,13]. The original FF of [8] is obtained if 7 = 0. The source 
variational principle is a one-parameter form that includes the potential energy and 
stress-displacement Reissner functionals as special cases [14-16]. 

(b) The ANDES variant of ANS is recovered if J 22 = .723 = 0 whereas j 2 i = a is a scaling 
parameter. The source variational principle is a one-parameter form that includes 
Reissner’s stress-displacement and Hu-Washizu’s functionals as special cases [18]. 

(c) If j 2 3 is nonzero, the last term in (2) may be viewed as being produced by a 
FF/ANDES combination. Such a combination remains unexplored. 




Box 2 Construction of the Basic Stiffness Matrix K& 


Step B.l . Assume a constant stress field, <r, inside the element. (This should be the 
element stress field that holds in the convergence limit; for structural elements the assumption 
would be on independent stress resultants. ) The associated boundary tractions are & n = 
<F.n, where n denotes the unit external normal on the element boundary 5. 

Step B.2. Assume boundary displacements, d, over S'. This field is described in terms of 
the visible element node displacements v (also called the connectors ) as 

d = v, (3) 

where Nj is an array of boundary shape functions. The boundary motions (3) must satisfy 
interelement continuity (or at least, zero mean discontinuity so that no energy is lost at 
interfaces) and contain rigid-body and constant-strain motions exactly. 

Step B.3. Construct the “lumping matrix” L that consistently “lumps” the boundary 
tractions <r n into element node forces, f, conjugate to v in the virtual work sense. That is, 

f = J N dn ff n dS = L<F. (4) 

In the above, N dn are boundary-system projections of conjugate to the surface tractions 

Step B.4' The basic stiffness matrix for a 3D element 

Kf, = tr 1 LEL t , (5) 

where E is the stress-strain constitutive matrix of elastic moduli, which are assumed to be 
constant over the element, and v = f v dV is the element volume measure. 

For a Kirchhoff plate bending element, stresses, strains and stress-strain moduli become 
bending moments, curvatures and moment-curvature moduli, respectively, and the integra- 
tion is performed over the element area A: 

K* = A~ l LDL t , (6) 

where D is the matrix of moment-curvature moduli. Specific examples for L are provided 
in Section 4. 
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Box 3 Construction of by the ANDES Formulation 


Step HA. Select reference lines (in 2D elements) or reference planes (in 3D elements) 
where “natural straingage” locations are to be chosen. By appropriate interpolation express 
the element natural strains e in terms of the “straingage readings” g at those locations: 

c=A c g, (7) 

where e is a strain field in natural coordinates that must include all constant strain states. 
(For structural elements the term “strain” is to be interpreted in a generalized sense.) 

Step H.2. Relate the Cartesian strains e to the natural strains: 


e = Te “ TA c g = Ag ' (8) 

at each point in the element. (If e = c, or if it is possible to work throughout in natural 
coordinates, this step is skipped.) 

Step H.3. Relate the natural straingage readings g to the visible degrees of freedom 

g-Qv, (9) 

where Q is a straingage- to- node displacement transformation matrix. Techniques for doing 
this vary from element to element and it is difficult to state rules that apply to every situation. 
.In the elements derived here Q is constructed by direct interpolation over the reference lines. 
(In general there is no unique internal displacement field u whose symmetric gradient is e 
or e, so this step cannot be done by simply integrating the strain field over the element and 
collocating u at the nodes.) 

Step H .4’ Split the Cartesian strain field into mean (volume-averaged) and deviatoric 

strains: 

e = e + e d = (A + A d )g, (10) 

where A = f v TA< dVjv, and ej = Ajg has mean zero value over V. This step may also 
be carried out on the natural strains if T is constant, as is the case for the elements here. 

Step H.5. The higher-order stiffness matrix is given by 

K>, = aQ T KjQ, with K d = f AjEAdrfT, (11) 

Jv 

where q ~ j 2 2 > 0 is a scaling coefficient (see Box 1). 

It is often convenient to combine the product of A and Q into a single strain-displacement 
matrix called (as usual) B, which splits into B and B<*: 

e =s AQv = (A + A^)Qv = (B + Bj)v = B v, (12) 

in which case 

K k = f BjEB d dV. (13) 

Jv 

The notation B c = A f Q is also used in the sequel. 



Box 4 Construction of K by the Conventional ANS Formulation 

Steps S.l to S.S. Identical to the first three steps H.l through H.S, in Box 3. The fourth 
step: strain splitting, is omitted. 

Step S.4 ■ The element stiffness matrix is given by 

K = Q T K a Q, with K a = [ A t EA dV. (14) 

Jv 

or, if B = AQ is readily available 


K = / B t EB dV. (15) 

Jv 

In general this stiffness matrix does not pass the individual element test of Bergan and 
Hanssen [6,7] (a strong form of the patch test that demands pairwise cancellation of node 
forces between adjacent elements in constant stress states). For this to happen, K must 
admit the decomposition 

K = K, + K h = v~'LEL t + K„, (16) 

where L is a force-lumping matrix derivable as per Box 2 and K /, is orthogonal to the rigid 
body and constant strain test motions. In other words, the ANS element must coalesce with 
the ANDES formulation with a = 1. The equivalence may be checked by requiring that 

B = AQ = u- 1 L t , (17) 

where A is the mean part of A. At the present there are no general techniques for explicit 
construction of strain fields that satisfy these conditions a priori. 

If the patch test is not satisfied, one should switch to the ANDES formulation by replacing 
the basic stiffness constructed from constant strain, namely uB EB, with one constructed 
from constant stress as in Box 2. 


2. THE TRIANGULAR ELEMENT 


2.1 Geometric Relations 


The geometry of an individual triangle is illustrated in Figure 1. The triangle has straight sides . Its geometry 
is completely defined by the location of its three corners, which are labelled 1,2,3, traversed counterclockwise. 
The element is referred to a local Cartesian system (x, y) which is usually taken with origin at the centroid 
0, whence the corner coordinates x>, yi satisfy the relations 

*i + *2 + *3 = 0, yi + V 2 + yz = 0- ( W 

Coordinate differences are abbreviated by writing x,j = x, — x ; -, and = t/,* — yj. The signed triangle area 
.4 is given by 


2.4 = 


1 

y\ 


l l 

*2 *3 

V2 2/3 


J 2ll/31 — *312/21 = *322/12 “ *122/32 = *132/23 “ *232/l3t 


(19) 


8 




Figure 1 The triangular element 


and we require that .4 > 0 . We shall make use of dimensionless triangular coordinates Ci» C 2 and ( 3 , linked 
by Cl + C 2 +C 3 = 1. The following well known relations between the triangular and Cartesian coordinates of 
a straight-sided triangle are noted for further use: 


x = x 1 C 1 + z 2 <2 + *3C3> y = yiCi + Z/2C2 + J/3C3, 


( 20 ) 


<» = 


1 

2 A 


XiVk - xtyj + (x - x 0 )yjk + (y - yo)xtj , 


(21) 


in which i, j and k denote positive cyclic permutations of 1, 2 and 3; for example, i = 2, j — 3, k — 1. (If 
the origin is taken at the centroid as in Figure 1, xq = j/o = 0.) It follows that 


2 A 
2 A 


<Ki 

dx 

dCi 

dy 


— 1/23, 


= *32, 


u £ = »'• 

o y 


2A ~d, X = yi21 
2A^ =*2i. 

dy 


( 22 ) 


Other intrinsic dimensions and ratios of use in future derivations are (see Figure 2) 


^ij — tji — \J x tj + Vij » c ij — j s ij — Vji/fiji 

a k — 2^4 1 ^ij , bij = XjJt *+* Vjiyki) / $-\j — * ^ij fyii (^-^) 

A ij = bij / iij = + yj%yki)/{ x ij H" y{j)i A ji = l — A,y = . 

Here = Cji is the length of side i-j and c, ; and s,y the cosine and sine, respectively, of angle (i — * 

Furthermore 6,y and 6j, are the projections of sides i-k and k-j } respectively, onto i-j ; A, ; and A ;> being 
the corresponding projection ratios. 

On each side i-j t define the dimensionless natur al c oordinates /i; ; as varying from 0 at i to 1 at j. The 
coordinate wj of a point not on the side is that of its projection on i-j. Obviously 


dx _ dy _ 

Wj ~ dua ~ Vji ' 


(24) 


9 



2.2 Displacements, Rotations, Curvatures 

As we are dealing with a Kirchhoff element, its displacement field is completely defined by the transverse 
displacement w(x } y) = u>(Ci> C 2 , <£ 3)1 positive upwards. In the present section we assume that w is unique 
and known inside the element; this assumption is relaxed later. The midplane (covariant) rotations about 
x and y are 9 X = dw/dy and 9 y = - dwjdx , respectively. Along side i-j with tangential direction t and 
external-normal n (see Figure 3) the tangential and normal rotations are defined as 


On = ^ = OxSij -0 y C ijy 
9% — ~~ — 0 £ Cij 4 * 9y Sij . 


(25) 


The visible degrees of freedom of the element collected in v (see Boxes 2-3) are 


V T = [uq 9 xl 9 y 1 W 2 9 x 2 0 y 2 W3 9 X 3 9 y 3 ]. (26) 

The first and second derivatives of the displacement w with respect to the Cartesian and triangular coordi- 
nates are linked by the relations (summation convention used) 


dw _ dw dC « _ 1 dw_ 

~di ~ iKiTi ~ 2AdO yjk ' 

dw _ dw d(i __ 1 dw 
dy dCi dy 2 A <90 


d 2 w d 2 w dQ dCj dw <9 2 0 _ 1 d 2 w 

~dx 2 ~ dUKj~fa dx + dCi dx 2 “ 4A 2 dQdCj yikyti ’ 
d 2 w _ d 2 w dQ d(j_ dw <9 2 Q 1 d 2 w 

dxdy dQdQ dx dy dQ dxdy AA 2 dCidCj^ k * * 

d 2 w d 2 w dQ dQ dw d 2 Q 1 d 2 w 

= dQdb'dy dCidCj XkjX ' k 


since d 2 (,?/dx 2 y d 2 ^ 2 /dxdy and d 2 (> 2 /dy 2 vanish on a straight-sided triangle, cf. Eq. (22). We can represent 
the second derivative relations in matrix form as 


d 2 w 

w 


y \ 3 *32 2x322/23 

T 

’ d 2 w ' 
d~ w 

1 

3/31 *13 2^132/31 

yjn *21 2x 2 li/l2 

< 

d 2 w 

’ “ 4A 2 

2 J/232/31 2X30X13 *322/31 + *132/23 


d 2 w 

; 


22/312/12 2x13X21 *132/12 4 >*212/31 

_ 2 y 122/23 2X21X32 *212/23 + *322/12 , 


w 1^2 

A 2 w 

d 2 w 

iKsOti 


or 


K — Ww^‘. 


( 30 ) 


The inverse relation does not exist. 
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2.3 Natural Curvatures 


The second derivatives of w with respect to the dimensionless side directions defined in Section 2.1 will be 
called the natural curvatures and denoted by \ij = d 2 w/d /i 2 j. Note that these curvatures have dimensions 
of displacement. The natural curvatures can be related to the Cartesian plate curvatures k xx = d 2 w/dx 2 , 
K yy = d 2 wfdy 2 and K xy = 2d 2 w/dxdy t by chain-rule application of (22): 



X = 


The inverse of this relation is 


Is 1 


> = 


T 1 

*21 

y \ i 

*211/21 


x 32 

Vl-i 

*321/32 

< 

*13 

y?3 

*13j/l3 „ 







) = 


1 

4 A 2 


31 


1/231/13 

*23*13 

J/23*31 + *321/13 


2/31 2/21 
*31*21 

2/31*12 + *13^21 


& 


^ 122/32 

*12*32 

J/l2*23 + *21 2/32 


= T" l ie. 


&Ml2 

®/ i 23 

d 2 w 

d/*3i 


(31) 


(32) 


or, in compact matrix notation 

* = T*. (33) 

A comparison of (29) with (31 )— (32) displays the advantages of natural curvatures over triangle-coordinate 
curvatures when the curvature field is to be constructed directly. On the other hand, (29) is useful when the 
transverse displacement w over the element is built as a function of the triangular coordinates. 


At this point we relax the requirement that the curvatures be derivable from a displacement field conse- 
quently the partial derivative notation will be discontinued. However, the foregoing transformations will be 
assumed to hold even if the curvature fields k and \ ar e not derivable from w. 


3. DIRECT CURVATURE INTERPOLATION 


3.1 The Straingage Readings 


ANS and ANDES plate bending elements are based on direct interpolation of natural curvatures. All elements 
discussed here adopt the three triangle sides as the reference lines defined in Box 3- The natural curvatures 
are assumed to vary linearly over each reference line, an assumption which is obviously consistent with cubic 
beam-like variations of it; over the sides. A linear variation on each side is determined by two straingage 
sample points, which we chose to be at the corners. 

Over each triangle side chose the isoparametric coordinates that vary from —1 at corner i to +1 at 
corner j. These are related to the fiij coordinates introduced in Section 2.1 by = 2 — 1. The Hermite 

interpolation of w over i-j is 


W= H(l-&>) 2 (2 + 6>) H-,(l (1 + e»'j ) 3 (2 — ft,) -U;(l+^) 2 (1 -£.,)] 


I 0 n) J 

where 9 n denotes the rotation about the external normal n on side ij. The natural curvature over side ij is 
given by 



Evaluating these relations at the nodes by setting = ±1 and converting normal rotations to x-y rotations 
through (25), we build the transformation 


' X12I1 ' 


'-6 

- 4 yai 

4 X 21 

6 

- 2 y 2 i 

2X21 

0 

0 

O' 

Xl2 1 2 


6 

2 y 2 i 

- 2 x 2 i 

-6 

4 y 2 i 

-4x 2 1 

0 

0 

0 

X23I2 


0 

0 

0 

-6 

— 4 y 3 2 

4x 32 

6 

- 2 y 32 

2 x 3 2 

X23I3 

► 

0 

0 

0 

6 

2 j /32 

— 2 X 32 

-6 

4 j /32 

-4x 32 

Yoi i ~ 


6 

— 2 yi 3 

2 xi 3 

0 

0 

0 

-6 

— 4 yi 3 

4x 13 

A 31 13 


.-6 

4 yi 3 

-4x 13 

0 

0 

0 

6 

2yi 3 

- 2 x 13 . 

• X 3 lli - 












r w i ' 

®x\ 

0y\ 

W2 

< #x2 > 

0y2 
W 3 

0x3 

< 0y3 > 


( 35 ) 


The left hand side is the natural straingage reading vector called g in Box 3 and thus we can express (35) as 


g = Qv. 


( 36 ) 


This relation holds for all elements discussed here. 

The six gage readings collected in g provide curvatures along the three triangle side directions at two corners. 
But nine values are needed to recover the complete curvature field over the element. The three additional 
values are the natural curvatures \23> \3i and X 12 at corners 1, 2 and 3, respectively. Three possibilities for 
the missing values are discussed below. 

3.2 The Average- Curvature Rule 

To each corner k assign the average natural curvature \ij °f the opposite side. This average is given by (34) 
evaluated at = 0. For example 

X 12 I 3 = X 12 I 1 + X 12 I 2 ) = 2/21 (0x2 — 0x1 ) + x 12(0y2 — 0y 1 )- (37) 

The natural curvature now can be interpolated linearly over the triangle: 


X12 = X12I1 Ci + X12I0 C2 + X12I3 C3 = X12I1 (Ci + 0C 3 ) + X 1 2 1 o (C2 + oC3)- ( 38 ) 

It is readily verified that under this rule the natural curvature \x 3 is constant over lines parallel to the 
triangle median that passes through node 3. Formulas for the other curvatures follow by cyclic permutation, 
from which we construct the matrix relation 




' X 12 1 

r Ci + 0 C 3 C 2 + 2 C 3 

0 0 

0 

0 


X23 ? = 

0 

0 

Ca + ki C 3 + 

0 C 1 

0 

0 

Ka. 

. X31 J 

0 

0 

0 c 

I 

C 3 + 0 C 2 

Cl + 


' 6 C 12 

( 3 C 21 — l)y2i 

(3Cl2 + 1)^21 

6 C 21 

(3C21 + 

i)y 2 i 

(3Cl2 

= 

0 

0 

0 

6C23 

(3C32 ~ 

i)y32 

(3C23 


. 6 C 13 

( 3 Ci 3 + l)yi 3 

(3C31 — 1)2/13 

0 

0 




0 

0 

0 






6 C 32 

(3C32 + l)y32 

(3C23 “ l)y32 

v, 





6 C 31 

(3ClS — l)yi 3 

(3C31 + l) x 13. 





= Ci “ C 2 , 

etc. In 

the notation of Box 3, 






( 39 ) 


X = A xo g = A X „Q v = B x 


m 
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where subscript a identifies the “averaging” rule (37), Since the natural curvatures vary linearly over the 
triangle, their mean values are obtained by evaluating (39) at the centroid £i = C2 = C3 — 1/3: 


[* 12 1 

'0 

-J/21 

*21 

0 

1/21 

— *21 

0 

0 

0 ■ 


*23 > = 

0 

0 

0 

0 

-1/32 

*32 

0 

1/32 

— *32 

V = 

l *31 J 

.0 

1/13 

-*13 

0 

0 

0 

0 

-yi3 

*13 . 



*x a 


(41) 


Finally, the Cartesian curvatures are given by 

k — TB xa v = B 0 v, (42) 

An explicit expression of these relations is easily obtained, but not required in what follows; however, that of 
the mean Cartesian curvatures 7c = TB xa v = B a v (a relation valid because T is constant over the triangle) 
is enlightening: 


] 

II 

K3| 

'0 

0 

2/32 

0 

0 

1/13 

0 

0 

y2i' 

1 *yy 

0 

Z32 

0 

0 

*13 

0 

0 

*21 

0 

[2K xy J 

0 

J/23 

*23 

0 

1/31 

*31 

0 

1/12 

*12 . 


v = B a v. 


( 43 ) 


3.3 The Projection Rule 

To each corner k assign the natural curvature \%j °f its projection onto the opposite side. This results in 
Xij being constant along lines normal to side ij. For equilateral triangles this agrees with the averaging 
rule, but not otherwise. The underlying motivation is to make the element insensitive to bad aspect ratios 
in cylindrical bending along side directions. 

To illustrate the application of this rule consider side 1-2. For node 3 take 

d^w 

X12I3 = Try = A12 xi2[i + A 2 i X12I2 1 ( 44 ) 

^12 3 

where A 12 and A 2 i are defined in (23). Proceeding similarly along the other sides we construct the matrix 
relation 


f \12 ’I 

Ci + ^ 12 C 3 

C 2 + -^21 Cs 

0 

0 

0 

0 



< X23 > = 

0 

0 

C 2 + ^23Cl 

C 3 + ^32Cl 

0 

0 

S’ 

(45) 

l X31 J 

0 

0 

0 

0 

£3 4- A 3 i£ 2 

Ci + A 13 C 2 _ 




or 


X = A xp g, k — TA XP g. 

where subscript p identifies the “projection” rule. As in the preceding rule, since T is constant we can do 
the strain-splitting step of Box 3 directly on the natural curvatures by evaluating at the centroid: 

A xp = (A xp 4- A x d p ) 

3(1 4" A12) 3(1 + A2O 0 0 0 0 

0 0 j(l 4- A03) 5(1 4- A 32 ) 0 0 

0 0 0 0 3(1 + ^31) §(1 + A 13 ) J 

C10 4 - A12C30 C20 + A21C30 0 0 0 0 

0 0 £20 4 - A03C10 C30 + A32C10 0 0 

0 0 0 0 (,30 4- A31C20 C10 4- A 13C20 J 


( 46 ) 


( 47 ) 


in 


which Qo = £* — Then 


Bp = TA XP Q = T(A X p 4 * Adp)Q = B p 4 - B<*p. 


( 48 ). 


The explicit expression of these matrices is not revealing and for the construction of the stiffness matrix 
presented in Appendix B it is better to leave (48) in product form. If all A coefficients are which happens 
for the equilateral triangle, the expressions reduce to those of the averaging rule. 
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3.4 The ‘Sliding Beam’ Rule 

This is a refinement of the average-curvature rule. Consider a fictitious beam parallel to side i — j sliding 
towards corner k. The end displacements and rotation of this beam are obtained by interpolating w cubically, 
9 n quadratically, and 9 t linearly, along sides i-k and j-k. Compute the mean natural curvature of this beam 
and assign to node it the limit as the beam reaches that corner. 

The required calculations can be simplified if we observe that the mean curvature of the sliding beam varies 
linearly as it moves from i-j, where it coincides with (41), to corner k. At one third of the way this mean is 
the natural centroidal curvature, which can then be readily extrapolated to k . These centroidal curvatures 
are given by x = B x ,v, where subscript s identifies the ‘sliding’ rule. A symbolic calculation yields the 
explicit form 


2 A 13 

” 2 (A 2 i 4 - A 31 ) 

2 A 12 

<* 2^13 

a 3 c 2 i 4 - a 2 Ci 3 

<* 3^21 

a 2 si 3 

a 3 s 2 i 4 - a 2 $i 3 

a 3 s 2 i 

2 A 23 

CM 

CM 

— 2 (Ai 2 4 - A 32 ) 

a l c 32 

a 3 c 2 l 

0^32 4 “ a 3 c 2 i 

a l 5 32 

^ 3 S 21 

a l 5 32 4 “ a 3 5 2 t 

- 2 (A 13 4 - A 23 ) 

2 A 31 

2 A 32 

fl 2 ci 3 4 - aic 32 

< 32^13 

<UC 32 

. a 2 st 3 4 - ais 32 

a 2 5 13 

ais 32 


where a f1 dj and are defined in Eqs. (23). Extrapolating to the opposite corners and interpolating over 
the triangle we get x = B xs v, with 


B 


T _ 
X * ~ 


6(— Ci + C2 + A 13C3) 
2t/2i(l — 3 Ci) + 3a 2 Ci 3 C 3 
2x2i(3Ci — 1) + 3a 2 si 3 C3 
6 (Ci - C2 4 * A03C3) 
2y 2 i(3C 2 “ 1) + 3aiC 3 2C3 
2x 2 i( 1 — 3£ 2 ) 4* 3ais 32 C3 
—6(A 2 3 + A 13 )C 3 
(3a 2 c 13 4- 3aiC 32 )C3 
(3a 2 si 3 4* 3ais 32 )C3 


— 6(A 2i 4- A 3 i)£i 
( 3a 3 c 2 i 4- 3a 2 ci 3 )Ci 
(3a 3 s 2 i + 3a 2 $i 3 )Ci 
6(— C 2 4- C 3 + A 2 iCi) 
22/32(1 — 3C 2 ) 4* 3a 3 c 2 iCi 
2x 32 (3^2 ”1)4" 3a 3 s 21 ^i 
6(Cs — C3 4- A 31 ^!) 
2y32(3C3 —1)4“ 3a 2 ci3^i 
2x 32 (l — 3£ 3 ) 4* 3a 2 si 3 £i 


6(C3 — Ci 4 * A 10 C 2 ) 
2 yi 3 ( 3 Ci — 1 ) 4 - 3 a 3 c 2 iC 2 
2xi 3 (l — 3Ci) 4“ 3a 3 s 2 iC 2 
— 6 ( A 1 2 4 - A 32 )C 2 
(3axc 32 4- 3a 3 c 2 i)C 2 
( 3 aiS 32 4 " 3 a 3 s 21 )C 2 
6 (— C 3 4 - Ci 4 - A 32 C 2 ) 
2 yi 3 (l ” 3 C 3 ) 4 - 3 a!C 32 C 2 
2xi 3 (3£ 3 ” 1) 4- 3ais 32 £ 2 . 


(50) 


It should be noted that Ax and Q are inextricably enmeshed in the above formula and cannot be easily 
separated. Premultiplication by T yields k ~ B,v. Evaluation of B, at the centroid yields B* = Lj /A, 
where = ATB XS is the force lumping matrix given in Eq. (56). 

A variation in the sliding-beam theme would consist of interpolating the normal rotation 9 n along i-k and j-k 
linearly rather than quadratically. This scheme turns out to be identical, however, to the average curvature 
rule and thus it provides nothing new. 

3.5 The Six Beam Lattice Rule 

In addition to the sides, consider three fictitious beams along the triangle medians. Determine the displace- 
ments and rotations at the triangle midpoints by the same interpolation procedure as in the sliding beam 
rule. The linear curvatures along the medians are thus readily computed. At each triangle corner we now 
know the curvatures in three directions: the two sides and the median. Consequently we can transform to 
x — y curvatures and interpolate these linearly over the element. This apparently new model gives, however, 
identical results to the projection rule, a result that can be a posteriori justified by geometric reasoning. 
Consequently this scheme will not be pursued further. 
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3.6 The ANS Elements 


Three ANS elements based on the previous interpolation rules may be constructed by following the pre- 
scription of Box 4. The stiffness matrices are identified as K a , K p , and K,, for averaging, projection, and 
sliding-beam, respectively. The following properties hold for these elements. 

Patch Test. Assuming that the element has constant thickness and material properties, K a and K, pass 
the individual element test, but K p does not. This claim can be analytically confirmed by applying the 
criterion of Eqs. (16)-(17), and noting that B a = L J /A and B, = L^/A, where L i and are the force 
lumping matrices derived in Section 4. 

Equivalence with DKT. K, turns out to be identical to the stiffness matrix of the Discrete Kirchhoff 
Triangle (DKT) element, which was originally constructed in a completely different way [3,4] that involves 
assumed rotation fields. Thus DKT is an ANS element , and also ( because of the equivalence noted below) an 
ANDES element . This equivalence provides the first variational justification of DKT, as well as the proof 
that DKT passes the patch test without any numerical verification. 

ANS/ANDES Equivalence. If the basic stiffness matrices K w and derived in Section 4.1 are used in 
conjunction with the averaging and sliding-beam rules, and a — 1, the ANDES formulation yields the same 
results as ANS if the element has constant thickness and material properties. (If the element has variable 
thickness, or the material properties vary, the equivalence does not hold.) The ANDES formulation used 
with the projection rule yields two elements, called ALR and AQR in the sequel, which differ in their basic 
stiffnesses. Both of these elements pass the patch test and are not equivalent to the ANS formulation. 


4. STIFFNESS MATRIX COMPUTATION 


4.1 The Basic Stiffness 


As explained on Box 2, the basic stiffness is obtained by constructing the lumping matrix L. In our case this 
is a 9 x 3 matrix that “lumps” an internal constant ben ding- moment field (7n xx , Tn yy , vn X y) to node forces 
f conjugate to v. 

On each element side, the constant moment field produces boundary moments rn nn and 7u n , referred to a 
local edge coordinate system n,t (see Figure 3): 



(51) 


The boundary motions d conjugate to m nn and m nt are dw/dn = —6 t and dw/dt = 0 n (see Figure 3). 
Given the degree of freedom configuration (25), the normal slope dwfdn = — 9 t along side i-j can at most 
vary linearly (it could be also taken as constant and equal to ^{9 tl + 9 t j) but the results are the same as for 
a linear variation). 

For the tangential slope (the rotation about the normal) dw/dt = 9 n there are three options: constant, 
linear and quadratic variation. But a constant 0 n = (wj — wi)/tij turns out to be equivalent to the quadratic 
variation and a constant 9 n = + #n>) equivalent to the linear variation. Consequently only the linear 

and quadratic cases need to be examined. 


Linear Normal Rotation. The variation of 9 t and 9 n along each side is linear: 




1-f 

0 


0 

i-i 


0 1+f 
0 0 


0 

l+£ 


Wi 

9u 

9m 

U 'j 

Oti 

@nj 


(52) 
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where £ = &j. Under this assumption one obtains [13] 


*0 0 yz2 0 0 yi3 0 0 V2\ 

L J = ^ 0 X 32 0 0 X 13 0 0 *21 0 , (53) 

|_0 2/23 *23 0 2/31 *31 0 2/12 * 12 . 

where superscript / stands for “linear The corresponding basic stiffness is 

K»j = A _1 LjDLjT, (54) 

where D is the Cartesian moment-curvature constitutive matrix resulting from the integration of E through 
the plate thickness. This matrix been used as component of the free formulation (FF) element presented in 
Ref. [13]. 

Quadratic Normal Rotation. A quadratic variation of 0 n can be accomodated in conjunction with the cubic 
variation of w along the side: 

' W 

9h 

r $ t \ _ i r o i-£ o o i-k o l J *ni 

\ 0 n ) ^ ~ 2 |.3(£ 2 ~ 1)/^ 0 + 1)(£ - 1) 3(^-1)/^ 0 i(3e-l)«+l)J] 

9 tj 
< &nj 

where £ = and £ = £ij. Then the resulting lumping matrix can be presented as 

— Cl2$i2 + c 31 5 31 — C31S31 -h Cl2^i2 ( 5 31 — c 3l) “ ( 5 12 ” C 12) 

2 ( s 12 x 12 + «31*3l) 2 (Ci 2 *12 + C 31 X 3l) cj 2 y 2 l + t /13 

— 2 ( S 12y21 + s 3ll/l3) — 2 ( c 1 2 1/21 + c 3lJ/l3) ~ 5 1 2 x 12 — s 31 X 31 

— C 2 3S 2 3 + Ci 2 Si 2 — Cl 2 Sl 2 + C 2 3S 2 3 (sj 2 CJ 2 ) _ ($ 2 3 — 0 23 ) 

2 ( s 1 2 x 12 + S 2 3 x 23) 2 ( c 12 x 12 + X 23 X 23 ) CJoJfri + c\ 3 X)32 (56) 

— §($122/21 + ^23^32) 2^(^122/21 + ^23^/23) ~ 5 12 X 12 - 5? 3 *23 

“ ^31 5 31 + ^ 23 ^23 “^23523 + ^31 5 31 (5§3 ~ ^23) ~ ( 5 31 “ C 3l) 

2( 5 23 x 23 + 5§! *3l) 2( C 23 X 23 + c 31 X 3l) ^3^32 + ^3l2/l3 

— ^( 5 23 2/32 + 53^13) — 2 (^23^32 + £31^13) ”-523X23 “ 531X31 

The corresponding basic stiffness matrix is denoted by 

K hq = A~ l L q DL[. (57) 

4.2 The Higher Order Stiffness 

The higher order stiffness for the ANDES elements described in Section 3 is 

K hx = aQ T K dx Q = qQ t A^DAj, d.d] Q = a J* Dj r DB Jx dA, (oS) 

where x = a,p, s for the average, projection and sliding-beam rules, respectively. (The last expression is 
appropriate when B & x is not easily factored into A^jQ, as in the sliding-beam rule.) Since varies 
linearly, if D is constant we could numerically integrate K <* x in (58) exactly with a three point Gauss rule; 
for example the three-midpoint formula. But as the element stiffness formation time is dominated b\ these 
calculations it is of interest to derive K /> in closed form. This is done in Appendix B for K hp , which from 
the numerical experiments appears to be the best performer. 
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5. NUMERICAL EXPERIMENTS: GENERAL DESCRIPTION 


An extensive set of numerical experiments has been run to assess the performance of the new ANDES 
elements based on the projection rule (ALR and AQR) and to compare them with other existing high- 
performance elements. Table 2 lists the tests, material properties and some relevant geometrical properties, 
whereas Table 3 lists elements, loading and mesh identifiers. 

An inspection of the element identifiers in Table 3 displays two important points: the difference in the results 
obtained with AQR and ALR can be attributed to their basic stiffness, whereas differences between AQR 
and DKT can be attributed to their higher order stiffness. With these facts in mind, we conducted first a 
set of distortion tests so that the less distortion sensitive combinations can be identified. Then, the best 
performers are submitted to a set of representative thin-plate bending problems in linear elasticity. 

The scaling a = 1.5 for ALR and a = 1.0 for AQR have been chosen to obtain energy balance in some simple 
cylindrical bending tests. No further adjustment of these parameters was made. In the distortion tests we 
included the results obtained with the free formulation (FF) element presented in [13], since that reference 
did not report such tests. 

Whenever the simply supported condition appears it implies that only the transverse displacement w is 
restrained. It is equivalent to the SSI condition described in Hughes’ textbook [21]. 

For tests involving an uniform distributed load q y two node-force computation schemes are usually reported: 

1. Triangular lumping (TL), in which one third of the load qA is assigned to each triangle corners, and 
nodal moments are set to zero. 


f T = 


qA 

3 


[1 0 0 1 0 0 1 0 0 ]. 


(59) 


2. Consistent lumping (CL), in which the element node force vector is 



£13 + Xl - 
^6 


2 1/23 + U22 ] . 


(60) 


This lumping was obtained using the transverse displacement w of the FF element in [13]. It is used for 
the ANS and ANDES elements as a matter of expediency, since for such elements a unique transverse 
displacement does not exist. 

Inasmuch as the present elements pass the linear patch test by virtue of their construction, no validation 
experiments along these lines are necessary once the elements are correctly programmed. 


6. DISTORTION TESTS 

6.1 Simply Supported Square Plate under Central Load 

This test was proposed by Kang [23]. The use of a coarse mesh exacerbates the distortion effect when far from 
of the converged solution. (In a fine mesh the distortion effect would be diluted.) The mesh and distortion 
parameter are shown In Figure 4. When the distortion parameter a approaches 2.5 the mesh converges to a 
four element cross- diagonal mesh. Results are reported as a percentage of the deterioration with respect to 
the undistorted mesh. 

The results given in Table 4 show that AQR is superior in this test. FF and ALR are the worst for (7 > 2. 
DKT and AQR display low deterioration rate from a = 2 up to a = 2.49, but DKT behaves poorly for a < 2. 
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Table 2. Key to Material and Geometrical Data 


Test 

Description 

Square plate 

Isotropic material v = 0, E = 1; thickness t = 1, plate span a = 10; 
load scaled so that center deflection w c = 100 

Cantilever beam 

Isotropic material v = 0, E = 1; thickness t = 1; load scaled so 
that center deflection w c = 100 

Twisted ribbon 

Isotropic material v = 0.25, E = 10 7 ; thickness t = 0.05; transverse 
load at tip so that Pb = — Pa = 1 

Rhombic cantilever 

Isotropic material v = 0.3, E = 10.5 10 6 ; thickness * = 0.125; 
uniform transverse load q = 0.26066 

Rhombic plate 

Isotropic material v — 0.3, E = 1; thickness t = 1, plate side 
a = 100, uniform transverse load q scaled so that w c — 100 


Table 3. Key to Element, Loading and Mesh Identifiers 


Key 

Explanation 


ANDES element Kw + 1.5K^ P 


ANDES element K^ -h K^ p 


FF element of [13] with 3-parameter scaling of K* 


ANS element Kt q + K*,: identical to DKT 

CL 

Consistent lumping (59) of uniform load q 

TL 

Triangular lumping (60) of uniform load q 

SDC 

In rhombic meshes, triangles obtained by splitting quadrilateral 
mesh units with short diagonal cuts 

LDC 

In rhombic meshes, triangles obtained by splitting quadrilateral 
mesh units with long diagonal cuts 
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6.2 Cantilever Beam 


A cantilever beam with a transverse load at the tip was selected for this test. Two meshes shown in Figure 5, 
A and B, are used to observe the effect of the element orientation under a linear bending state. The results 
are reported in Table 5. Also shown in this table is the ratio of the computed tip deflection to the exact 
value w ex for zero distortion. 

For mesh A, AQR is the best performer closely followed by DKT. FF and ALR behave poorly. 

For mesh B FF is the best performer in terms of deterioration, followed by AQR, DKT and ALR. However 
it must be noted that FF and ALR recover only 77% of the exact solution. This is a serious drawback in 
elements supposedly capable of providing an appropriate response for linear bending. This shortcoming can 
be attributed to the basic stiffness K*/ which is the same for both elements. AQR and DKT recover almost 
99% of the response for both meshes. 

6.3 Twisted Ribbon 

This test has been selected to assess the distortion effect under a field which combines bending and twisting. 
The test uses mesh B of Figure 5. The results shown in Table 6 indicate that AQR and DKT are the least 
distortion sensitive elements for this problem. 


20 



undistorted 


distorted 


Figure 4. Square plate: mesh for distortion analysis 


Table 4* Distortion Analysis: Centrally Loaded SS Square Plate: Percent 
Error of Center Deflection with Respect to Undistorted Mesh 


Element 


Distortion parameter 


type 

0.50 

1.00 

1.50 

2.00 

2.49 

ALR 

0.83 

2.65 

5.05 

7.88 

10.38 

AQR 

0.17 

-0.14 

-1.59 

-3.29 

-4.40 

DKT 

-0.95 

-3.46 

-6.29 

-8.06 

-8.42 

FF 

0.81 

2.27 

3.69 

4.85 

-13.50 
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Table 5. Distortion Analysis: Cantilever Beam 
Percent Error at Node C with Respect to Undistorted Mesh 


Mesh 

Element 

Distortion parameter 

wc/wctx 


type 

1.00 

3.00 

4.90 

(no distortion) 

A 

ALR 

— 10.70 

-19.80 

8.40 

1.031 

A 

AQR 

0.15 

0.10 

-2.05 

0.995 

A 

DKT 

0.20 

-0.59 

3.41 

0.982 

A 

FF 

-7.75 

-17.30 

-18.35 

0.974 

B 

ALR 

0.20 

3.00 

45.90 

0.764 

B 

AQR 

-0.10 

0.40 

-2.85 

0.995 

B 

DKT 

-0.13 

-1.09 

-3.49 

0.979 

B 

FF 

-0.05 

-0.15 

2.20 

0.769 
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Table 6. Distortion Analysis - Twisted Ribbon 
Loss of Symmetry under Distortion (Mesh B) 


Element 

type 

Node 

Distortion parameter 
1.00 3.00 4.90 

ALR 

A 

1.016 

1.122 

1.363 


B 

1.013 

1.098 

1.076 

AQR 

A 

0.989 

0.966 

0.945 


B 

1.010 

1.029 

0.995 

DKT 

A 

0.993 

0.978 

0.940 


B 

1.006 

1.015 

1.018 

FF 

A 

0.983 

0.933 

0.789 


B 

0.994 

0.877 

0.877 


7 . CONVERGENCE STUDIES 

From the distortion test results, it can be concluded that elements whose basic stiffness is are less 
distortion sensitive. Consequently only results for the AQR and DKT elements are presented in the following 
studies. 

7*1 Square Plate 

In this analysis a square plate with either simply-supported or clamped edges is considered. Due to symmetry 
only one quarter of the plate is modeled. The two different mesh orientations, A and B, used in the analysis 
are illustrated in Figure 6. The number of elements used is 2 N 2 y where N is the number of side subdivisions. 

For the cases involving a concentrated load, Figures 7 and 8 show that for both meshes AQR converges faster 
and is less sensitive to mesh orientation than DKT. 

In the case of uniform loading with triangular lumping, Figures 9 and 10, the convergence is uniform for 
all the meshes and elements. For the simply-supported condition all answers are within the 5% error limit 
for 1 V = 4. Clearly DKT converges faster in this case. For the clamped condition and N = 4, DKT(A) is 
outside the 5% error limit. 

For consistent force lumping, the results shown in Table 7 indicate a dramatic improvement of AQR. DKT 
also improves in the sense that becomes less mesh sensitive and that all the results are within 5% error for 
N = 4. 

7*2 Rhombic Cantilever 

The test involves a rhombic cantilevered plate subjected to uniform load. This problem was presented in [3] 
to test the DKT element with reference given to an experimental deflection result; however, no convergence 
analysis was performed. This has been done here taking into account the two possible mesh subdivision 
patterns, SDC and LDC, depicted in Figure 11. Triangular force lumping has been used. 

The results are shown in Table 8. For the LDC mesh DKT converges from above to an answer 4% below 
the experimental value quoted in [3]. On the other hand, AQR converges from below. For the SDC mesh 
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Mash A 


Mash B 


C 


L 


Figure 6. Meshes for square plate convergence studies 


both elements behave identically and converge to a value 4% under the experimental one. 

It is clear from these results that the experimental tip deflection given in [3] is in error by about -H % with 
respect to the analytical value for the material properties quoted. The apparently small error for the 2x2 
DKT/LDC mesh is thus fortuitous. 

7.3 Simply Supported Rhombic Plate 

This problem poses severe difficulties for ordinary finite element methods because of the presence of a 
singularity in the bending moments at the obtuse corner. A detailed description of this problem may be seen 
for example in [13]. The acute angle a = 30° was selected for the test. Again both SDC and LDC meshes 
were tried. 

The results are shown in Table 9. For the SDC meshes AQR and DKT show slight difference and almost 
the same rate of convergence. For the LDC meshes DKT is too flexible whereas AQR converges faster. 



Figure 7. Central deflection of centrally loaded SS square plate 



Figure 8. Central deflection of centrally loaded clamped square plate 
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Figure 9. Central deflection of uniformly loaded SS square plate 
with TL force lumping 



Figure 10. Central deflection of uniformly loaded clamped square plate 
with TL force lumping 
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Table 7. Uniformly Loaded Square Plate with CL Force Lumping: 
Percent Error of Central Deflection 


Support 

Element 

type 

Mesh 

type 

Mesh over quarter plate 
lxl 2x2 4x4 8x8 

ss 

DKT 

A 

31.73 

4.49 

1.01 

0.24 



B 

4.55 

5.37 

1.56 

0.41 


AQR 

A 

16.28 

2.20 

0.47 

0.11 



B 

-L55 

2.30 

0.74 

0.20 

Clamped 

DKT 

A 

46.35 

14.90 

4.10 

1.03 



B 

-21.60 

2.08 

1.30 

0.36 


AQR 

A 

26.65 

8.26 

1.87 

0.44 



B 

-41.20 

-3.22 

-0.28 

-0.05 


Table 8. Rhombic Cantilever: Percent Difference 
of Tip- A Deflection with Respect to Experimental Value of [3] 


Mesh 

Element 

Subdivision of whole plate 

type 

type 

4x4 

8x8 

16 x 16 

LDC 

DKT 

2.3 

-3.7 

-4.0 


AQR 

-17.8 

— 10,4 

-6.0 

SDC 

DKT 

-6.7 

-5.0 

-4.0 


AQR 

-6.3 

-5.0 
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Figure 11. Rhombic cantilever: meshes for convergence studies 


Table 0. Uniformly Loaded SS Rhombic Plate with TL Force 
Lumping: Percent Error in Center Deflection 


Mesh 

Element 

Subdivision of 

whole plate 

type 

type 

4x4 

8x8 

16 x 16 

SDC 

DKT 

11.05 

4.07 

2.86 


AQR 

13.86 

4.56 

2.89 

LDC 

DKT 

80.97 

22.64 

7.51 


AQR 

6.85 

-0.36 

-2.91 
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8. CONCLUSIONS 


The main conclusions of the present study can be summarized as follows. 

1. The ANDES formulation represents a variant of the ANS formulation that merits serious study. The 
key advantages of ANDES over ANS are: 

(a) a priori satisfaction of the patch test. Although this advantage is less clear for elements where 
ANS and ANDES coalesce for constant thickness and material properties, it reappears for more 
general cases. 

(b) The separation of the higher order stiffness allows the application of a scaling parameter. Fur- 
thermore it opens the possibility for the energy-balanced combination with other formulations as 
per Eq. (2), although this possibility presently remains unexplored. 

2. The study of plate bending elements shows that the widely used DKT element is both an ANS and 
ANDES element. This discovery provides a variational foundation hereto lacking and analytically proves 
(because of the ANDES connection) that DKT passes the patch test. 

3. The numerical results clearly demonstrate that the choice of basic stiffness is of paramount importance in 
the behavior of elements based on the ANDES formulation. Of the two elements sharing the quadratic- 
rotation basic stiffness, namely AQR and DKT, the former has excelled in geometric distortion tests 
and in convergence studies that involve concentrated forces. For other cases the performance of AQR 
and DKT is similar, and superior to those elements that use the linear-rotation basic stiffness. 

The numerical experiments have not addressed questions of material sensitivity such as element perfor- 
mance for highly anisotropic and composite plates. This behavior, as well as the possibility of applying this 
technology to C° bending elements, is currently under investigation. 
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Appendix A. SANITIZING INCOMPATIBLE ELEMENTS 


The stiffness-splitting technique summarized in Box 1 provides a systematic way for “sanitizing” existing 
nonconforming bending elements that do not pass the patch test. The technique amounts to the replacement 
of the basic stiffness. The main steps will be briefly outlined for the simplest such element: the BCIZ triangle 
proposed in 1965 by Bazeley et al [5]. The assumed transverse displacement is given explicitly in [12] as 


w = < 
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where t/j = yu — J/31 1 5/2 = 1/23 - 1/12, 5/3 = J/31 — J/23. *1 — *21 — *13, *2 = *32 — *21. *3 — *13 — *32- The 
strain-displacement matrix B is obtained by double differentiation with respect to the triangular coordinates 
and application of (30): 
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in which W is given by (29), and 
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Split the strain-displacement equations as 

K = K + Kj = (B + Bd)v, 
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where B = Bq 4 i(Bi + Bj 4 B 3 ) , B<* = B — B. Then the “sanitized” stiffness matrix is 


K = K b + a j bJD 6 B d dA, (65) 

where K& is one of the basic stiffness matrices derived in Section 4.1. The free formulation leads to the 
same result but in a less direct manner, because w would have to be decomposed into rigid body, constant 
curvature and higher order states. Although the corrected element passes the patch test it is unlikely to be 
competitive with ANDES elements in distortion insensitivity as this property appears to depend on relaxing 
curvature compatibility conditions. 
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Appendix B. EXPLICIT REPRESENTATION OF HIGHER ORDER STIFFNESS 


To obtain an explicit representation of K^ P) begin by defining 

C n C12 C13 
C22 C23 
C23 


C = T DT = 


(66) 

[symm 

which can be interpreted as a constitutive matrix that relates the natural moments T T m to the natural 
curvatures x * Then 
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Carrying out the congruential transformation K* p = Q T K < i p Q with MACSYMA yields 
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The same stiffness expression applies for K^ a , if one sets Ai 2 = A 23 = A 3i = 
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